TD3 - 17/05/21 - Simulation et Monte Carlo¶
- Exercice 7 : MCMC
- Exercice 9 : Cross-Entropy Method
Exercice 7 : MCMC¶
Rappel :
7.1.¶
- Quelle loi reconnaissez-vous ?
7.2.¶
- Calculer les lois conditionnelles de chaque composante.
- Mettre en oeuvre le Gibbs sampler correspondant
Code
In [3]:
from scipy import stats
import numpy as np
def rand_gauss_gibbs(n, rho):
# initialisation : nous pouvons choisir autre valeur
X = np.zeros((n, 2))
# itération gibbs sampler
for t in range(n-1):
X[t, 0] = stats.norm.rvs(rho*X[t, 1], np.sqrt(1-rho**2), 1)
X[t+1, 1] = stats.norm.rvs(rho*X[t, 0], np.sqrt(1-rho**2), 1)
# ne pas oublier la dernière simulation de X1
X[n-1, 0] = stats.norm.rvs(rho*X[n-1, 1], np.sqrt(1-rho**2), 1)
return(X)
rn = rand_gauss_gibbs(10, 0.6)
In [4]:
rn
Out[4]:
- Représenter dans le plan : des contours de la densité cible et l'évolution de la chaîne de Markov simulée
In [5]:
import matplotlib.pyplot as plt
# fonction qui affiche les contours de la densité cible et l'échantillon simulé
def plot_gibbs_contour(n, rho):
# simuler X1 et X2
rn = rand_gauss_gibbs(n, rho)
# contours de la densité cible
rv = stats.multivariate_normal(np.zeros(2), np.array([1, rho, rho, 1]).reshape(2, 2))
x, y = np.mgrid[min(rn[:, 0]):max(rn[:, 0]):.1, min(rn[:, 1]):max(rn[:, 1]):.1]
pos = np.dstack((x, y))
plt.contour(x, y, rv.pdf(pos))
plt.title(f"rho = {rho}")
plt.xlabel(f"X1")
plt.ylabel(f"X2")
# représentation de la chaine de Markov simulée
plt.scatter(rn[:,0], rn[:,1], color = "steelblue")
return None
In [6]:
plot_gibbs_contour(n = 500, rho = 0.5)
- Que se passe-t-il quand $\rho \to 1$ ?
In [7]:
n = 100
rhos = [0, 0.2, 0.5, 0.7, 0.9, 0.99]
plt.figure(figsize = (16, 10))
for i, rho in enumerate(rhos):
plt.subplot(2, 3, i+1)
plot_gibbs_contour(n = n, rho = rho)
lorsque $\rho \to 1$, $X_2 = X_1$ p.s.
rq : on peut montrer (et voir) facilement que lorsque $\rho \to -1$, $X_2 = -X_1$ p.s.
Nous utilisons les package python statsmodels pour évaluer la performance MCMC
In [8]:
from statsmodels.graphics.tsaplots import plot_acf
import pandas as pd
n = 5000
def plot_gibbs_ACF(n, rho, dim = 0):
# simuler X1 et X2
rn = rand_gauss_gibbs(n, rho)
# plot
plot_acf(rn[:,dim])
plt.title(f"Autocorrelation pour X{dim+1}")
return None
def plot_gibbs_trace(n, rho, dim = 0):
# simuler X1 et X2
rn = rand_gauss_gibbs(n, rho)
# plot
plt.plot(np.arange(n), rn[:,dim])
plt.title(f"Trace X{dim+1}")
return None
In [9]:
rho = 0.5
plot_gibbs_ACF(n, rho, dim = 0)
plot_gibbs_ACF(n, rho, dim = 1)
plt.show()
In [10]:
rho = 0.5
plt.figure(figsize=(16, 4))
plt.subplot(1,2,1)
plot_gibbs_trace(n, rho, dim = 0)
plt.subplot(1,2,2)
plot_gibbs_trace(n, rho, dim = 1)
plt.show()
Que se passe-t-il si $\rho$ proche de 1 ?
In [11]:
plot_gibbs_ACF(n, 0.99, dim = 0)
plot_gibbs_ACF(n, 0.99, dim = 1)
plt.show()
In [12]:
rho = 0.99
plt.figure(figsize=(16, 4))
plt.subplot(1,2,1)
plot_gibbs_trace(n, rho, dim = 0)
plt.subplot(1,2,2)
plot_gibbs_trace(n, rho, dim = 1)
plt.show()
7.3.¶
- Programmer un algorithme RWHM (Hastings-Metropolis avec marche aléatoire)
In [13]:
def RWMH(n, rho, sigma):
X = np.zeros((n,2))
for t in range(n-1):
# proposition
x = stats.multivariate_normal(X[t,:], sigma).rvs(1)
rvn = stats.multivariate_normal(np.zeros(2), np.array([1, rho, rho, 1]).reshape(2, 2))
if np.log(np.random.rand()) < np.log(rvn.pdf(x)) - np.log(rvn.pdf(X[t,:])):
X[t+1,:] = x
else:
X[t+1,:] = X[t,:]
return X
rho = 0.8
tau = 0.8
X = RWMH(n, rho, tau * np.eye(2))
X[:10]
Out[13]:
In [37]:
# fonction qui affiche les contours de la densité cible et l'échantillon simulé
def plot_RWMH_contour(n, rho, tau, sigma = np.eye(2)):
# simuler X1 et X2
rn = RWMH(n, rho, tau*sigma)
# contours de la densité cible
rv = stats.multivariate_normal(np.zeros(2), np.array([1, rho, rho, 1]).reshape(2, 2))
x, y = np.mgrid[min(rn[:, 0]):max(rn[:, 0]):.1, min(rn[:, 1]):max(rn[:, 1]):.1]
pos = np.dstack((x, y))
plt.contour(x, y, rv.pdf(pos))
plt.title(f"rho = {rho}, tau = {tau}")
plt.xlabel(f"X1")
plt.ylabel(f"X2")
# représentation de la chaine de Markov simulée
plt.scatter(rn[:,0], rn[:,1], color = "steelblue")
return None
def plot_RWMH_ACF(n, rho, tau, sigma = np.eye(2), dim = 0):
# simuler X1 et X2
rn = RWMH(n, rho, tau*sigma)
# plot
plot_acf(rn[:,dim])
plt.title(f"Autocorrelation pour X{dim+1} avec rho = {rho}, tau = {tau}")
return None
def plot_RWMH_trace(n, rho, tau, sigma = np.eye(2), dim = 0):
# simuler X1 et X2
rn = RWMH(n, rho, tau*sigma)
# plot
plt.plot(np.arange(len(rn)), rn[:,dim])
plt.title(f"Trace X{dim+1} avec rho = {rho}, tau = {tau}")
return None
- voici ci-dessous les paramètres que nous allons analyser :
In [15]:
taus = [0.001, 1, 5, 100]
rhos = [0.2, 0.9, 0.99]
In [16]:
plt.figure(figsize=(16, 4))
for i, tau in enumerate(taus):
plt.subplot(1, 4, i+1)
plot_RWMH_contour(n = n, rho = rhos[0], tau = tau)
In [17]:
plt.figure(figsize=(16, 4))
for i, tau in enumerate(taus):
plt.subplot(1, 4, i+1)
plot_RWMH_contour(n = n, rho = rhos[1], tau = tau)
In [18]:
plt.figure(figsize=(16, 4))
for i, tau in enumerate(taus):
plt.subplot(1, 4, i+1)
plot_RWMH_contour(n = n, rho = rhos[2], tau = tau)
In [19]:
plot_RWMH_ACF(n, rhos[2], taus[1], sigma = np.eye(2), dim = 0)
plot_RWMH_ACF(n, rhos[2], taus[1], sigma = np.eye(2), dim = 1)
plt.show()
In [20]:
plt.figure(figsize=(16, 4))
plt.subplot(1,2,1)
plot_RWMH_trace(n, rhos[2], taus[1], sigma = np.eye(2), dim = 0)
plt.subplot(1,2,2)
plot_RWMH_trace(n, rhos[2], taus[1], sigma = np.eye(2), dim = 1)
plt.show()
Ne tenons plus compte des doublons :
In [34]:
def RWMH(n, rho, sigma):
X = [[0, 0]]
for t in range(n-1):
# proposition
x = stats.multivariate_normal(np.array(X[-1]), sigma).rvs(1)
rvn = stats.multivariate_normal(np.zeros(2), np.array([1, rho, rho, 1]).reshape(2, 2))
if np.log(np.random.rand()) < np.log(rvn.pdf(x)) - np.log(rvn.pdf(X[-1])):
X.append(x)
return np.array(X)
In [39]:
plot_RWMH_ACF(n, rhos[2], taus[1], sigma = np.eye(2), dim = 0)
plot_RWMH_ACF(n, rhos[2], taus[1], sigma = np.eye(2), dim = 1)
plt.show()
In [41]:
plt.figure(figsize=(16, 4))
plt.subplot(1,2,1)
plot_RWMH_trace(n, rhos[2], taus[1], sigma = np.eye(2), dim = 0)
plt.subplot(1,2,2)
plot_RWMH_trace(n, rhos[2], taus[1], sigma = np.eye(2), dim = 1)
plt.show()
Bilan des courses :
Gibbs sampling est un bon outil pour construire des chaînes de Markov
Cependant lorsque les v.a marginales sont très corrélées ($\rho \to 1$) nous avons des anomalies : loi simulée différente de la loi cible
nous pouvons améliorer la performance de la simulation en utilisant un algorithme RWHM (où $\Sigma$ de la loi auxilière $\mathcal{N}(\mu, \Sigma)$ doit être bien choisie, pas "trop grand", pas "trop petit").
Exercice 9 : Cross-Entropy¶
Code :
In [ ]:
from numpy.random import multivariate_normal
def rosenbrock(x):
y = 0
for i in range(len(x)-1):
y = y + 100 * (x[i+1]-x[i])**2 + (x[i]-1)**2
return(y)
def optCE(fun = rosenbrock, n = 1000, d = 3, eps = 1e-8, max_iter = 1000, tau = 0.01):
"""
this function compute the arg max of the function 'fun' by cross-entropy method
"""
mu = np.random.rand(d)
sigma = 10*np.eye(d)
t = 0
while True:
t += 1
# step 1 : sample Yi
Y = multivariate_normal(mean = mu, cov = sigma, size=n)
# step 2 : pick the top Yi that minimize the function
samples = - np.array(list(map(fun, Y)))
n_top = round(tau * n)
top_ind = np.argsort(samples)[::-1][:n_top]
Y_top = Y[top_ind, :]
# step 3 : estimate by MLE mu and sigma
mu = np.mean(Y_top, 0)
sigma = np.cov(Y_top.T)
# step 4 : stopping criter
if np.max(sigma) < eps or t >= max_iter:
print(t)
print(mu)
break
return mu
In [ ]:
optCE(n = 1000, d = 3)
In [ ]:
optCE(n = 1000, d = 5)
In [ ]:
optCE(n = 1000, d = 10)
In [ ]:
optCE(n = 10000, d = 10)
Bilan des courses :
- la méthode CE peut-être utilisé soit pour une estimation soit pour une optimisation
pour une optimisation, on doit tenir en compte plusieurs paramètres :
- si $n$ (taille de l'échantillon) devient grand alors l'optimisation est plus précise
- si $d$ (dimension) devient grand alors l'optimisation est plus compliqué
l'initialisation de $\theta_0$ est bien sûr très important, ici $\theta_0 = (\mu_0, \sigma^2_0)$ et :
- idéalement $\mu_0$ doit être proche de la vraie valeur (sinon on choisit une valeur déterministe ou aléatoire)
- $\sigma_0$ doit être assez grand afin de considérer toutes les possibilités